clear 
est clear

********************************************************************************
*	Programs
********************************************************************************

* program for the OLS regressions... 
capture program drop lm_regs
program define lm_regs, eclass

	*** Syntax 
	syntax varlist(numeric) [, y_bl(varlist numeric) ctrls(varlist numeric) treatment(varlist numeric) rounds(numlist) t_suff(string)]

	*** Unadjusted Regressions 
qui{
	foreach y of local varlist {
		foreach r in `rounds' {
		*** Create the centered block variables 
			capture drop blocks*
			qui tab block, gen(blocks)
			qui foreach v of varlist blocks* {
					replace `v' = . 			if missing(`y')
					su `v' 						if round == `r' & !missing(`y')
					replace `v' = `v' - r(mean) if round == `r' & !missing(`y')
				}
		*** Linear Model: treatment, Block FE, Cluster SEs at Club level, weights (tracking) 
		eststo `y'_`t_suff'_`r':		reg `y' b0.`treatment'##c.(blocks*) [pw=w] if round==`r', vce(cluster club)
			qui {
				*** add scalars
				* control mean
				qui su `y' [aw=w] if `treatment'==0 & round==`r'   
				estadd scalar c_mean = r(mean) 	
				estadd scalar c_sd = r(sd)
				* p-value 
				test 1.`treatment' == 0 
				estadd scalar pval = r(p) 
			}
		}
	} 
}
	**** Adjusted Regressions 
qui {
	foreach y of local varlist {
		foreach r in `rounds' {
		*** Create the centered block variables 
			capture drop blocks*
			qui tab block, gen(blocks)
			qui foreach v of varlist blocks* {
					su `v' 						if round == `r' & !missing(`y')
					replace `v' = `v' - r(mean) if round == `r' & !missing(`y')
				}
		*** Linear Model: treatment, Block FE, Cluster SEs at Club level, weights (tracking), BL covars fully interacted
		eststo `y'_`t_suff'_`r'_adj:	reg `y' b0.`treatment'##c.(`ctrls' `y_bl' blocks*) [pw=w] if round==`r', vce(cluster club)
			qui {
				*** add scalars
				* control mean
				qui su `y' [aw=w] if `treatment'==0 & round==`r'   
				estadd scalar c_mean = r(mean) 	
				estadd scalar c_sd = r(sd)
				* p-value 
				test 1.`treatment' == 0 
				estadd scalar pval = r(p) 
			}
		}
	}
}
end



* program to extract q values for the outcomes specified BY round
capture program drop get_qvals
program define get_qvals, eclass

	*** Syntax 
	syntax varlist(numeric) [, rounds(numlist) adj(string) t_suff(string) method(string) iv(string)]

	*** if LATE change the suffix 
	if ("`iv'"=="yes") local t_suff = "iv"
	di "`t_suff'"
	
	*** Create vector of p-values
	* matrix to store them
	local varcount : word count `varlist'
	matrix pvals_tab = J(`varcount',1,.) 
	matrix colnames pvals_tab = pvals
	* loop over outcomes and rounds
	qui foreach r in `rounds' {
	    local i = 0
		foreach y of local varlist {
		    local ++i
			* restore estimates
			est restore `y'_`t_suff'_`r'`adj' 
			est replay
			* collect pvals in matrix 
			if ("`iv'"=="yes") {
				local beta = e(b)[1,"1.treat_ther"]
				local beta_se = sqrt(e(V)["1.treat_ther","1.treat_ther"])
				local beta_z = `beta'/`beta_se'
				local beta_p = 2*(1-normal(abs(`beta_z')))
				matrix pvals_tab[`i',1] = `beta_p'
				matrix list pvals_tab
			}
			else if ("`iv'"!="yes") {
				matrix list r(table)
				matrix pvals_tab[`i',1] = r(table)["pvalue","1.treat_`t_suff'"]
			}
			* save in a dataset 
			capture drop pvals`r'
			svmat pvals_tab, names(col)
			rename pvals pvals`r'
			
			*** get q-values  
			capture drop qvals`r'
			qqvalue pvals`r', method(`method') qvalue(qvals`r') 
		}
	}
	*** add to estimates 
	* loop over outcomes and rounds
	qui foreach r in `rounds' {
	    local i = 0
		foreach y of local varlist {
		    local ++i
			* restore estimates
			est restore `y'_`t_suff'_`r'`adj' 
			* add as a scalar
			local qv = qvals`r' in `i'
			di `qv'
			estadd scalar qval = `qv'
		}
	} 
end











********************************************************************************
*	I - Impacts of Therapy vs. Control
********************************************************************************

use "${data}/Analysis_wide.dta", clear
distinct block club

*** is it easier in long mode? I think it might be...
keep se_score* ed_enrolled* competencies* everpregnant* newlypregnant* evermarried* newlymarried* riskysex* desired_fert* time_pref* paid_work* life_exp* ed_aspiration* wantpreg* wantmarr* age_b everpregnant_b evermarried_b PPI_score_b phq8_score_b w? block cr_id club treat treat_ther
reshape long ed_enrolled competencies se_score everpregnant newlypregnant evermarried newlymarried riskysex desired_fert time_pref paid_work life_exp ed_aspiration wantpreg wantmarr w, i(block cr_id club treat ) j(round)

*** drop +cash folks in longer follow ups 
drop if treat==2 & round>=2
tab round 

*** Run regs 
* a bunch we can do normally 
foreach y in ed_enrolled newlypregnant newlymarried riskysex time_pref paid_work desired_fert {
	lm_regs `y', rounds(1 2 3) ///
		treatment(treat_ther) t_suff("ther") y_bl(phq8_score_b) ///
		ctrls(age_b everpregnant_b evermarried_b PPI_score_b) 
}
* competencies missing at 12 months 
lm_regs competencies ed_aspiration wantpreg wantmarr,  rounds(1 3) ///
	treatment(treat_ther) t_suff("ther") y_bl(phq8_score_b) ///
	ctrls(age_b everpregnant_b evermarried_b PPI_score_b)
* self efficacy not available at baseline 
lm_regs se_score,  rounds(2 3) ///
	treatment(treat_ther) t_suff("ther") y_bl(phq8_score_b) ///
	ctrls(age_b everpregnant_b evermarried_b PPI_score_b)
* life expectancy at RR only
lm_regs life_exp,  rounds(1) ///
	treatment(treat_ther) t_suff("ther") y_bl(phq8_score_b) ///
	ctrls(age_b everpregnant_b evermarried_b PPI_score_b)

	
*** compute q-values 
* "simes" is the Benjamini Hochberg (1995) step-up procedure
* primary outcomes 
get_qvals ed_enrolled competencies 			newlypregnant newlymarried riskysex, rounds(1) adj(" ") t_suff("ther") method("simes")
get_qvals ed_enrolled competencies 			newlypregnant newlymarried riskysex, rounds(1) adj("_adj") t_suff("ther") method("simes")
get_qvals ed_enrolled 			   se_score newlypregnant newlymarried riskysex, rounds(2) adj(" ") t_suff("ther") method("simes")
get_qvals ed_enrolled 			   se_score newlypregnant newlymarried riskysex, rounds(2) adj("_adj") t_suff("ther") method("simes")
get_qvals ed_enrolled competencies se_score newlypregnant newlymarried riskysex, rounds(3) adj(" ") t_suff("ther") method("simes")
get_qvals ed_enrolled competencies se_score newlypregnant newlymarried riskysex, rounds(3) adj("_adj") t_suff("ther") method("simes")
* secondary outcomes 
get_qvals time_pref paid_work desired_fert, rounds(1 2 3) adj(" ") t_suff("ther") method("simes")
get_qvals time_pref paid_work desired_fert, rounds(1 2 3) adj("_adj") t_suff("ther") method("simes")


** empty regression for competencies in round 2? 
eststo competencies_empty: reg cr_id block

*** unadjusted table
* Midline
#delimit ;
esttab 	ed_enrolled_ther_2 competencies_empty se_score_ther_2 newlypregnant_ther_2 newlymarried_ther_2 riskysex_ther_2 
		time_pref_ther_2 paid_work_ther_2 desired_fert_ther_2
	using "${tables}/results-hc-both2.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_ther) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-value")) 
	; 
#delimit cr 
* Endline
#delimit ;
esttab 	ed_enrolled_ther_3 competencies_ther_3 se_score_ther_3 newlypregnant_ther_3 newlymarried_ther_3 riskysex_ther_3 
		time_pref_ther_3 paid_work_ther_3 desired_fert_ther_3
	using "${tables}/results-hc-both3.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_ther) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-value")) 
	; 
#delimit cr 

*** adjusted table
* Midline
#delimit ;
esttab 	ed_enrolled_ther_2_adj competencies_empty se_score_ther_2_adj newlypregnant_ther_2_adj newlymarried_ther_2_adj riskysex_ther_2_adj 
		time_pref_ther_2_adj paid_work_ther_2_adj desired_fert_ther_2_adj
	using "${tables}/results-hc-both2-adj.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_ther) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-value")) 
	; 
#delimit cr 
* Endline
#delimit ;
esttab 	ed_enrolled_ther_3_adj competencies_ther_3_adj se_score_ther_3_adj newlypregnant_ther_3_adj newlymarried_ther_3_adj riskysex_ther_3_adj 
		time_pref_ther_3_adj paid_work_ther_3_adj desired_fert_ther_3_adj
	using "${tables}/results-hc-both3-adj.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_ther) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-value")) 
	; 
#delimit cr 




*** DO the R1 only outcomes table 
est clear

* regs 
* life expectancy at RR only
lm_regs competencies time_pref ed_aspiration paid_work desired_fert wantpreg wantmarr life_exp,  rounds(1) ///
	treatment(treat_ther) t_suff("ther") y_bl(phq8_score_b) ///
	ctrls(age_b everpregnant_b evermarried_b PPI_score_b)


* RR secondary outcomes 
get_qvals competencies time_pref desired_fert paid_work , rounds(1) adj(" ") t_suff("ther") method("simes")
get_qvals competencies time_pref desired_fert paid_work , rounds(1) adj("_adj") t_suff("ther") method("simes")
get_qvals ed_aspiration wantpreg wantmarr life_exp, rounds(1) adj(" ") t_suff("ther") method("simes")
get_qvals ed_aspiration wantpreg wantmarr life_exp, rounds(1) adj("_adj") t_suff("ther") method("simes")

 
* RR Only secondary outcomes 
* unadjusted
#delimit ;
esttab 	desired_fert_ther_1
		time_pref_ther_1  
		paid_work_ther_1 
		competencies_ther_1
		ed_aspiration_ther_1
		wantpreg_ther_1 
		wantmarr_ther_1 
		life_exp_ther_1
	using "${tables}/results-hc-sec-ther-r1only.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_ther) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-value")) 
	; 
#delimit cr 


* adjusted 
#delimit ;
esttab 	desired_fert_ther_1_adj
		time_pref_ther_1_adj  
		paid_work_ther_1_adj 
		competencies_ther_1_adj
		ed_aspiration_ther_1_adj
		wantpreg_ther_1_adj 
		wantmarr_ther_1_adj 
		life_exp_ther_1_adj
	using "${tables}/results-hc-sec-ther-r1only_adj.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_ther) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-value")) 
	; 
#delimit cr 

















**********************************************************************
*	Phase I - Impacts of Cash + Therapy over just Therapy
**********************************************************************


use "${data}/Analysis_wide.dta", clear
distinct block club

*** is it easier in long mode? I think it might be...
keep se_score* ed_enrolled* competencies* newlypregnant* newlymarried* riskysex* desired_fert* time_pref* paid_work* life_exp* ed_aspiration* wantpreg* wantmarr* age_b everpregnant_b evermarried_b PPI_score_b phq8_score_b w? block cr_id club treat treat_ther treat_cash
reshape long ed_enrolled competencies se_score newlypregnant newlymarried riskysex desired_fert time_pref paid_work life_exp ed_aspiration wantpreg wantmarr w, i(block cr_id club treat ) j(round)

*** drop controls folks in longer follow ups 
drop if treat==0 
drop if round<=1
tab round 

*** value label of cash arm 
label define treat_cash 0 "IPT-G" 1 "IPT-G+", modify
 
*** Run regs 
* a bunch we can do normally 
foreach y in ed_enrolled se_score newlypregnant newlymarried riskysex time_pref paid_work desired_fert {
	lm_regs `y', rounds(2 3) ///
		treatment(treat_cash) t_suff("cash") y_bl(phq8_score_b) ///
		ctrls(age_b everpregnant_b evermarried_b PPI_score_b) 
}
* competencies missing at 12 months 
lm_regs competencies,  rounds(3) ///
	treatment(treat_cash) t_suff("cash") y_bl(phq8_score_b) ///
	ctrls(age_b everpregnant_b evermarried_b PPI_score_b)

*** compute q-values 
* "simes" is the Benjamini Hochberg (1995) step-up procedure
* primary outcomes 
get_qvals ed_enrolled 			   se_score newlypregnant newlymarried riskysex, rounds(2) adj(" ") t_suff("cash") method("simes")
get_qvals ed_enrolled 			   se_score newlypregnant newlymarried riskysex, rounds(2) adj("_adj") t_suff("cash") method("simes")
get_qvals ed_enrolled competencies se_score newlypregnant newlymarried riskysex, rounds(3) adj(" ") t_suff("cash") method("simes")
get_qvals ed_enrolled competencies se_score newlypregnant newlymarried riskysex, rounds(3) adj("_adj") t_suff("cash") method("simes")
* secondary outcomes 
get_qvals time_pref paid_work desired_fert, rounds(2 3) adj(" ") t_suff("cash") method("simes")
get_qvals time_pref paid_work desired_fert, rounds(2 3) adj("_adj") t_suff("cash") method("simes")

** empty regression for competencies in round 2? 
eststo competencies_empty: reg cr_id block

*** unadjusted table
* Midline
#delimit ;
esttab 	ed_enrolled_cash_2 competencies_empty se_score_cash_2 newlypregnant_cash_2 newlymarried_cash_2 riskysex_cash_2 
		time_pref_cash_2 paid_work_cash_2 desired_fert_cash_2
	using "${tables}/results-hc-both2-cash.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_cash) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("IPT-G mean" "IPT-G SD" "Observations" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
* Endline
#delimit ;
esttab 	ed_enrolled_cash_3 competencies_cash_3 se_score_cash_3 newlypregnant_cash_3 newlymarried_cash_3 riskysex_cash_3 
		time_pref_cash_3 paid_work_cash_3 desired_fert_cash_3
	using "${tables}/results-hc-both3-cash.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_cash) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("IPT-G mean" "IPT-G SD" "Observations" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 

*** adjusted table
* Midline
#delimit ;
esttab 	ed_enrolled_cash_2_adj competencies_empty se_score_cash_2_adj newlypregnant_cash_2_adj newlymarried_cash_2_adj riskysex_cash_2_adj 
		time_pref_cash_2_adj paid_work_cash_2_adj desired_fert_cash_2_adj
	using "${tables}/results-hc-both2-cash-adj.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_cash) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("IPT-G mean" "IPT-G SD" "Observations" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
* Endline
#delimit ;
esttab 	ed_enrolled_cash_3_adj competencies_cash_3_adj se_score_cash_3_adj newlypregnant_cash_3_adj newlymarried_cash_3_adj riskysex_cash_3_adj 
		time_pref_cash_3_adj paid_work_cash_3_adj desired_fert_cash_3_adj
	using "${tables}/results-hc-both3-cash-adj.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_cash) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("IPT-G mean" "IPT-G SD" "Observations" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
